
. /* ------------------------------------------------------------------------ ** 
>     C. Wunder, A. Wiencierz, J. Schwarze, H. K�chenhoff:    
>     Well-being over the life span: semiparametric evidence from British 
>     and German longitudinal data
>     
>     ESTIMATION OF SEMIPARAMETRIC REGRESSIONS
> 
>     Data organization:  person-year observations
>     Data source:        German Socio-Economic Panel Study (SOEP), years 1986-2007.
> ** ------------------------------------------------------------------------ */ 
. 
. /*  --------( load data )-------------------------------------------------- */
. version 10

. matrix  drop _all

. use     "SOEP.dta", clear
(PanelWhiz v2.0 Nov 2007 <john@panelwhiz.eu>)

. rename  age x

. rename  gls_life_sat gls_y

. 
. /*  --------( set knots )-------------------------------------------------- */
. mkspline    spline 16 = x, marginal displayknots pctile

             |     knot1      knot2      knot3      knot4      knot5      knot6      knot7 
-------------+-----------------------------------------------------------------------------
           x |        23         27         30         33         36         39         42 

             |     knot8      knot9     knot10     knot11     knot12     knot13     knot14 
-------------+-----------------------------------------------------------------------------
           x |        45         48         52         55         59         63         68 

             |    knot15 
-------------+-----------
           x |        74 

. global      N_knots = r(N_knots)            

. global      N_spl = r(N_knots)+1            

. matrix      knots = r(knots)'

. 
. /*  --------( define macros & prepare data )------------------------------ */
. global  xvar ///
>         female      disabled    hospital    educ        ln_netto    ln_hhsize ///
>         german      fulltime    parttime    unempl      single      divorced  ///
>         widowed     west        d_year4         attr_in_1   attr_in_2   attr_in_3 ///
>         d_year5     d_year6     d_year8     d_year9     d_year11    d_year12  ///
>         d_year13    d_year14    d_year15    d_year16    d_year17    d_year18  ///
>         d_year19    d_year20    d_year21    d_year22    d_year23    d_year24  ///
>         one

. 
. tokenize ${xvar}

. 
. global  gls_xvar ///
>         gls_`1'  gls_`2'  gls_`3'  gls_`4'  gls_`5'  gls_`6'  gls_`7'  gls_`8'  ///
>         gls_`9'  gls_`10' gls_`11' gls_`12' gls_`13' gls_`14' gls_`15' gls_`16' ///
>         gls_`17' gls_`18' gls_`19' gls_`20' gls_`21' gls_`22' gls_`23' gls_`24' ///
>         gls_`25' gls_`26' gls_`27' gls_`28' gls_`29' gls_`30' gls_`31' gls_`32' ///
>         gls_`33' gls_`34' gls_`35' gls_`36' gls_`37'

. 
. global copy_gls_xvar ///
>         _gls_`1'  _gls_`2'  _gls_`3'  _gls_`4'  _gls_`5'  _gls_`6'  _gls_`7'  _gls_`8'  /
> //
>         _gls_`9'  _gls_`10' _gls_`11' _gls_`12' _gls_`13' _gls_`14' _gls_`15' _gls_`16' /
> //
>         _gls_`17' _gls_`18' _gls_`19' _gls_`20' _gls_`21' _gls_`22' _gls_`23' _gls_`24' /
> //
>         _gls_`25' _gls_`26' _gls_`27' _gls_`28' _gls_`29' _gls_`30' _gls_`31' _gls_`32' /
> //
>         _gls_`33' _gls_`34' _gls_`35' _gls_`36' _gls_`37'

. 
. global  TPF ///
>         spline2     spline3     spline4     spline5     spline6     ///
>         spline7     spline8     spline9     spline10    spline11    ///
>         spline12    spline13    spline14    spline15    spline16    

. foreach var of global TPF {
  2.     qui replace `var' = `var'^3
  3. }

. gen spline1_sq = spline1^2

. gen spline1_cu = spline1^3

.                                
. global spline_vars spline1 spline1_sq spline1_cu ${TPF} 

. foreach var of global spline_vars{
  2.     qui bysort persnr: egen im_`var' = mean(`var')
  3.     qui gen gls_`var' = `var' - theta_i*im_`var'
  4. }

. drop im_*

. foreach var of global spline_vars{
  2.     qui gen _gls_`var' = gls_`var'
  3. }

. 
. /*  --------( estimation of semiparametric model eq. in 4 )----------------- */
. xtmixed gls_y gls_spline1 gls_spline1_sq gls_spline1_cu ${gls_xvar}, noconst /// 
>         || _all: gls_spline2-gls_spline$N_spl, noconstant cov(identity) ///
>         , technique(bfgs) emdots

Performing EM optimization: ....................

Performing gradient-based optimization: 

Iteration 0:   log restricted-likelihood =  -425054.6  
Iteration 1:   log restricted-likelihood =  -425054.6  

Computing standard errors:

Mixed-effects REML regression                   Number of obs      =    253044
Group variable: _all                            Number of groups   =         1

                                                Obs per group: min =    253044
                                                               avg =  253044.0
                                                               max =    253044


                                                Wald chi2(40)      = 384871.65
Log restricted-likelihood =  -425054.6          Prob > chi2        =    0.0000

--------------------------------------------------------------------------------
         gls_y |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
---------------+----------------------------------------------------------------
   gls_spline1 |  -.6724533   .4228479    -1.59   0.112     -1.50122    .1563133
gls_spline1_sq |   .0268438     .01839     1.46   0.144       -.0092    .0628876
gls_spline1_cu |  -.0003679   .0002665    -1.38   0.167    -.0008901    .0001544
    gls_female |   .0741481   .0145509     5.10   0.000     .0456287    .1026674
  gls_disabled |  -.4519398   .0143441   -31.51   0.000    -.4800537   -.4238259
  gls_hospital |  -.0120578    .000328   -36.76   0.000    -.0127007   -.0114149
      gls_educ |   .0338148   .0024974    13.54   0.000     .0289199    .0387097
  gls_ln_netto |   .4919393   .0095259    51.64   0.000     .4732688    .5106098
 gls_ln_hhsize |  -.1942077   .0119516   -16.25   0.000    -.2176324   -.1707831
    gls_german |   .0529215    .020209     2.62   0.009     .0133126    .0925304
  gls_fulltime |   .0789905   .0111111     7.11   0.000     .0572132    .1007678
  gls_parttime |   .0191999   .0124057     1.55   0.122    -.0051149    .0435146
    gls_unempl |  -.5970624    .013975   -42.72   0.000    -.6244528    -.569672
    gls_single |  -.1737692   .0165677   -10.49   0.000    -.2062412   -.1412972
  gls_divorced |   -.137162   .0183242    -7.49   0.000    -.1730767   -.1012473
   gls_widowed |  -.1960326   .0228048    -8.60   0.000    -.2407291   -.1513361
      gls_west |   .5105448   .0165043    30.93   0.000      .478197    .5428927
   gls_d_year4 |  -.1731851   .0194437    -8.91   0.000     -.211294   -.1350762
 gls_attr_in_1 |  -.1085155   .0131398    -8.26   0.000     -.134269    -.082762
 gls_attr_in_2 |  -.1255516    .013125    -9.57   0.000    -.1512761    -.099827
 gls_attr_in_3 |  -.0828003   .0117849    -7.03   0.000    -.1058982   -.0597023
   gls_d_year5 |   -.241719   .0197453   -12.24   0.000    -.2804192   -.2030189
   gls_d_year6 |  -.2441144   .0200025   -12.20   0.000    -.2833186   -.2049101
   gls_d_year8 |  -.0517673   .0204122    -2.54   0.011    -.0917744   -.0117602
   gls_d_year9 |  -.2700781   .0190799   -14.16   0.000     -.307474   -.2326822
  gls_d_year11 |  -.3986364   .0194745   -20.47   0.000    -.4368056   -.3604671
  gls_d_year12 |  -.3863249    .019743   -19.57   0.000    -.4250204   -.3476294
  gls_d_year13 |  -.3935779     .01979   -19.89   0.000    -.4323657   -.3547902
  gls_d_year14 |  -.5152072   .0197962   -26.03   0.000     -.554007   -.4764074
  gls_d_year15 |  -.4379885   .0200943   -21.80   0.000    -.4773726   -.3986044
  gls_d_year16 |  -.4019374   .0203307   -19.77   0.000    -.4417849     -.36209
  gls_d_year17 |  -.4588068   .0202451   -22.66   0.000    -.4984866   -.4191271
  gls_d_year18 |  -.4321291    .020554   -21.02   0.000    -.4724143    -.391844
  gls_d_year19 |   -.228007   .0190213   -11.99   0.000    -.2652879    -.190726
  gls_d_year20 |  -.2829232   .0192237   -14.72   0.000    -.3206008   -.2452455
  gls_d_year21 |  -.4338474    .019251   -22.54   0.000    -.4715787   -.3961161
  gls_d_year22 |  -.3126028   .0194733   -16.05   0.000    -.3507697    -.274436
  gls_d_year23 |  -.4179277   .0196757   -21.24   0.000    -.4564914    -.379364
  gls_d_year24 |  -.3928076   .0198675   -19.77   0.000    -.4317473   -.3538679
       gls_one |   8.813081   3.236981     2.72   0.006     2.468715    15.15745
--------------------------------------------------------------------------------

------------------------------------------------------------------------------
  Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]
-----------------------------+------------------------------------------------
_all: Identity               |
   sd(gls_s~e2..gls_s~16)(1) |   .0002839   .0000879      .0001547    .0005209
-----------------------------+------------------------------------------------
                sd(Residual) |    1.29718   .0018236      1.293611    1.300759
------------------------------------------------------------------------------
LR test vs. linear regression: chibar2(01) =   555.58 Prob >= chibar2 = 0.0000
(1) gls_spline2 gls_spline3 gls_spline4 gls_spline5 gls_spline6 gls_spline7 gls_spline8
    gls_spline9 gls_spline10 gls_spline11 gls_spline12 gls_spline13 gls_spline14
    gls_spline15 gls_spline16

. est store model1b

. global N_fe = e(k_f)

. 
. /*  --------( prediction: smooth function )-------------------------------- */
. predict reff*, reffects level(_all)     

. global all_vars ${xvar} ${spline_vars}

. foreach var of global all_vars{     
  2.     qui replace gls_`var' = `var'   
  3. }

. replace gls_one = 1
(0 real changes made)

.                        
. local i = 1                      

. local j = `i' + 1

. while `i' < $N_spl {
  2. gen re_product`i' = reff`i'*spline`j'
  3. local i = `i' + 1
  4. local j = `i' + 1
  5. }

. egen blups_rcoef = rowtotal(re_product1-re_product$N_knots)

. drop re_product*

. adjust ${gls_xvar}, generate(yfixed)    

-------------------------------------------------------------------------------------------
     Dependent variable: gls_y     Equation: gls_y     Command: xtmixed
       Created variable: yfixed
   Variables left as is: gls_spline1, gls_spline~q, gls_spline~u
 Covariates set to mean: gls_female = .51646749, gls_disabled = .10785871,
                         gls_hospital = 1.8507256, gls_educ = 11.545087,
                         gls_ln_netto = 7.9742905, gls_ln_hhs~e = .93653946,
                         gls_german = .83492989, gls_fulltime = .44324307,
                         gls_parttime = .1218049, gls_unempl = .06990089,
                         gls_single = .19605286, gls_divorced = .06776687,
                         gls_widowed = .06582254, gls_west = .77702692,
                         gls_d_year4 = .03595817, gls_attr_i~1 = .06747443,
                         gls_attr_i~2 = .10576026, gls_attr_i~3 = .13411502,
                         gls_d_year5 = .03468962, gls_d_year6 = .03373326,
                         gls_d_year8 = .03258722, gls_d_year9 = .04685351,
                         gls_d_year11 = .04508702, gls_d_year12 = .04373152,
                         gls_d_year13 = .04483015, gls_d_year14 = .04635953,
                         gls_d_year15 = .04415042, gls_d_year16 = .04320592,
                         gls_d_year17 = .04704715, gls_d_year18 = .04518582,
                         gls_d_year19 = .07169109, gls_d_year20 = .06972305,
                         gls_d_year21 = .07437047, gls_d_year22 = .0715014,
                         gls_d_year23 = .06802769, gls_d_year24 = .06511121, gls_one = 1
-------------------------------------------------------------------------------------------

----------------------
      All |         xb
----------+-----------
          |   -4.77291
----------------------
     Key:  xb  =  Linear Prediction

. gen EBLUP = blups_rcoef + yfixed

. 
. /*  --------( 95% Bonferroni-corrected CIs )------------------------------- */
. sort x

. matrix accum CTC = _gls_spline1 _gls_spline1_sq _gls_spline1_cu     /// 
>     ${copy_gls_xvar} _gls_spline2-_gls_spline$N_spl, noconstant 
(obs=253044)

. global N_fe_spl = $N_fe + $N_knots      

. matrix nullcol = J($N_knots,$N_fe,0)

. matrix nullrow = J($N_fe,$N_fe_spl,0)

. matrix D = (nullrow \ nullcol,I($N_knots))

. matrix estimates = e(b)             

. matrix est_sig_e = estimates[1,"lnsig_e:_cons"]

. matrix est_sig_u = estimates[1,"lns1_1_1:_cons"]

. scalar sigma_e = exp(est_sig_e[1,1])

. scalar sigma_u = exp(est_sig_u[1,1])

. scalar lambda = (sigma_e^2/sigma_u^2) 

. matrix Cov = sigma_e^2*inv(CTC + lambda*D)  

. collapse (mean) EBLUP spline1 spline1_sq spline1_cu spline2-spline$N_spl ${xvar}, by(x)

. mkmat spline1 spline1_sq spline1_cu ${xvar} spline2-spline$N_spl, matrix(C_x)

. matrix Var = C_x*Cov*C_x'

. matrix Var_f = vecdiag(Var)

. matrix Var_f = Var_f'

. svmat Var_f, names(Var_f)

. gen se_f = sqrt(Var_f)

. scalar z = invnormal(1-0.025/83)

. gen ci_upper = EBLUP + z*se_f

. gen ci_lower = EBLUP - z*se_f

. 
. /*  --------( plot of smooth function )------------------------------------ */
. cap svmat knots, names(knots)           

. levelsof knots, local(myknots)
23 27 30 33 36 39 42 45 48 52 55 59 63 68 74

. twoway (rarea ci_upper ci_lower x, bcolor(gs12)) ||     ///
>      (line EBLUP x, sort connect(L) lstyle(solid)) ,    ///
>      legend(off) yscale(range(3 8)) ylabel(3(1)8)       ///
>      ytitle(life satisfaction) xtick(20(5)100) ytick(3(0.5)8)   ///
>      xmtick(`myknots', tpos(inside) tlength(*2) tlc(black))

. qui log close
